home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byt0486a.arc / ATOMTLBX.LBR / ATSPGM.FOR < prev    next >
Text File  |  1986-04-11  |  4KB  |  150 lines
  1.  
  2. c*+*+*+*+*+*
  3. c     This program was produced by the    ATOMCC    translator version 7.10
  4. c                    Copyright (C) 1985, Y. F. Chang
  5. c*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+
  6. c Portions (c) Copyright, Microsoft Corp., 1981.  All rights reserved.
  7. c This program was written for the following inputs
  8. c
  9. C FIRST PAINLEVE TRANSCENDENT
  10. C     DIFF(Y,T,2) = 6.0*Y*Y + T
  11. c--------
  12. c no instructions in second input block
  13. c--------
  14.       COMMON /IPASS/ LENSER,LENVAR,MPRINT,MSTIFF,LRUN,
  15.      + KTRDCV,KNTSTP,KTSTIF,KXPNUM,KDIGS,KENDFG,NTERMS,NOPT
  16.      A /RPASS/ RADIUS,ERRLIM,ADJSTF,RCREAL,RCIMAG,RDCERR
  17.      B /CPASS/ START,END,ORDER
  18.      C /DPASS/ H,HNEW,XPRINT,DLTXPT
  19.       DIMENSION TMPS(  36,  1)
  20.       CHARACTER*6 NAMES
  21.       EQUIVALENCE (TMPS(1,1),Y(1))
  22.       DIMENSION NAMES(1), Y(36), T(2), TMPAAB(30), TMPAAA(30)
  23.       DATA NAMES(1)/'Y.....'/
  24.       Y(33) = 1.1
  25.    10 FORMAT(72H   ATOMCC  Ver. 7.10, Copyright (C) 1985, Y. F. Chang; S
  26.      Aolution results./9H   ******)
  27.    11 FORMAT(/5X,11HStep number,I6,13H at the point,1P1E12.4/1X,
  28.      A 9Hvalues of )
  29.    12 FORMAT(1X, A6,(1X,1P4E13. 5))
  30.    13 FORMAT(5X,21HStepsize adjusted to ,1PE13.5)
  31.    14 FORMAT(/5X,35HThe solution stopped normally after, I4,24H steps as
  32.      a set by nsteps. )
  33.    16 FORMAT(/5X,63HThe adjustment for stepsize seems to be in a loop. P
  34.      Alease try a /5X,22Hshorter series length. )
  35.       WRITE(*,10)
  36. c--------
  37. c Initialize variables to default values.
  38. c--------
  39.       NSTEPS = 40
  40.       H = 1.E0
  41.       ERRLIM = 1.E- 6
  42.       LENSER = 30
  43.       MPRINT = 4
  44.       NTERMS = 2
  45.       KTRDCV =    1
  46.       ADJSTF = 1.E-2
  47.       MSTIFF =    0
  48.       DLTXPT = 0.E0
  49. c--------
  50. c start of third  input block
  51. c--------
  52. C READ INTEGRATION INTERVAL AND INITIAL CONDITIONS.
  53.       READ(5,1010) START,END,Y(1),Y(2)
  54.  1010 FORMAT(4F10.3)
  55.       WRITE(*,1020) START,END,Y(1),Y(2)
  56.  1020 FORMAT(' SOLVE THE FIRST PAINLEVE TRANSCENDENT' /
  57.      + ' INTERVAL:          ',2F10.3 /
  58.      + ' INITIAL CONDITIONS:',2F10.3 /)
  59. c--------
  60. c end of third    input block
  61. c--------
  62. c More initializations
  63. c--------
  64.       DLTXPT = SIGN(DLTXPT,(END-START))
  65.       H = SIGN(H,(END-START))
  66.       KDIGS =  6
  67.       XPRINT = START + DLTXPT
  68.       KXPNUM = 35
  69.       LENVAR = 36
  70.       LRUN = 1
  71.       KTSTIF = 0
  72.       NUMEQS =     1
  73.       IF(LENSER.GT.(LENVAR- 6)) LENSER = LENVAR -  6
  74.       IF(MPRINT.LT.2) GO TO 17
  75.       WRITE(*,11) KTSTIF,START
  76.      K = Y(33)
  77.      WRITE(*,12) NAMES(K),Y(1), Y(2)
  78. c--------
  79. c Loop for integration steps.  Inside the loop, print the desired output
  80. c--------
  81.    17 DO 27 KINTS=1,NSTEPS
  82.      KOUNT = 0
  83.      KNTSTP = KINTS
  84.    19     CONTINUE
  85.      T(1) = START
  86.      T(2) = H
  87.      Y(2) = Y(2)*(H)
  88. c--------
  89. c Preliminary series calculations
  90. c--------
  91.      TMPAAA(1) = 6.E0*Y(1)
  92.      TMPAAB(1) = TMPAAA(1)*Y(1)
  93.      Y(3) = (TMPAAB(1) + T(1))*(H*H/2.E0)
  94.      TMPAAA(2) = 6.E0*Y(2)
  95.      TMPAAB(2) = TMPAAA(1)*Y(2) + TMPAAA(2)*Y(1)
  96.      Y(4) = (TMPAAB(2) + T(2))*(H*H/6.E0)
  97. c--------
  98. c Loop for series calculations
  99. c--------
  100.      DO 23 K= 5,LENSER
  101.         KA = K - 1
  102.         KB = K - 2
  103.         TMPAAA(KB) = 6.E0*Y(KB)
  104.         TMPAAB(KB) = 0.E0
  105.         KZ = 1 + KB
  106.         DO 30 N=1, KB
  107.            L = KZ - N
  108.            TMPAAB(KB) = TMPAAB(KB) + TMPAAA(N)*Y(L)
  109.    30         CONTINUE
  110.         Y(K) = (TMPAAB(KB))*(H*H/(KB*KA))
  111. c--------
  112. c Test and adjust H to avoid over/under flow.
  113. c--------
  114.         IF(MSTIFF.GE.20 .AND. KTSTIF.GT.0) GO TO 23
  115.         TMP = ABS(Y(K))
  116.         IF(TMP.LE.1.E-35) GO TO 23
  117.         IF(TMP.LT.1.E20 .AND. TMP.GT.1.E-20) GO TO 23
  118.         IF(KTSTIF.NE.0 .AND. TMP.LT.1.0) GO TO 23
  119.         KOUNT = KOUNT + 1
  120.         IF(KOUNT.LT.9) GO TO 22
  121.         WRITE(*,16)
  122.         GO TO 28
  123.    22        CONTINUE
  124.         Y(2) = Y(2)/(H)
  125.         H = H * TMP**(0.3/(1-K))
  126.         IF(MPRINT.GE.4) WRITE(*,13) H
  127.         GO TO 19
  128.    23     LRUN = 1
  129. c--------
  130. c Calculate radius of convergence and take optimum step.
  131. c--------
  132.      CALL RDCV(TMPS,LENVAR,NUMEQS,NAMES)
  133.    24     CALL RSET(TMPS,LENVAR,NUMEQS,NAMES)
  134. c--------
  135. c no instructions in fourth input block
  136. c--------
  137.    25     GO TO (26,28,24), KENDFG
  138.    26     H = SIGN(RADIUS,H)
  139.      START = START + HNEW
  140.      IF(MPRINT.LT.4) GO TO 27
  141.      WRITE(*,11) KNTSTP, START
  142.      K = Y(33)
  143.      WRITE(*,12) NAMES(K),Y(1), Y(2)
  144.    27 CONTINUE
  145.       WRITE(*,14) NSTEPS
  146.    28 CONTINUE
  147.    29 STOP
  148.       END
  149. K = Y(33)
  150.      WRITE(*,12) NAMES(K),Y(1),